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We review some computer algorithms for the simulation of off-lattice clusters 
grown from a seed, with emphasis on the diffusion-limited aggregation, ballistic ag- 
gregation and Eden models. Only those methods which can be immediately extended 
to distinct off-lattice aggregation processes are discussed. The computer efficiencies 
of the distinct algorithms are compared. 
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INTRODUCTION 



Growth processes occurring far from equilibrium are widespread in nature and technoL 



neurite formation ^\. Computer models for the growth of clusters, generally constituted 
of identical particles, are useful tools for the understanding of aggregation phenomena. 
The main contribution of such models is to provide pathways to investigate the underly- 
ing physical ingredients ruling the properties observed in growth phenomena. One of the 
most intriguing features of the fractal structures found in nature and computer models is 
the scale invariance emerging without fine-tuning of any parameter, in contrast with usual 



ogy. Examples include electrodeposition viscous fingering [2], bacterial colonies [3t], and 



critical phenomena in which scale invariance only emerges at a critical point 



a. 



The foremost example of nonequilibrim growt 
(DLA) model introduced by Witten and Sander 



1 model is the diffusion-limited aggregation 



in 1981. The rules of the DLA model are 
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based on an iterative stochastic process in which the particles, one at a time, follow Brownian 
trajectories until they touch and stick in an aggregate. Despite its simple rules, the DLA 
model leads to very complex aggregates with multiscale properties [3, Q] and multifractality 
in the growth-site probability distribution jol. llol|. 

If the random walks in the DLA model are replaced by ballistic trajectories at random 
directions, we have the ballistic aggregation (BA) model [l2| proposed by Void to describe 
colloid aggregation. Differently from DLA, the BA model generates asymptotically non- 
fractal clusters (fractal dimension equal to the space dimension) characterized by a power 
law approach to the asymptotic regime Q, [l^ . 

Finally, a third standard aggregation process was proposed by Eden 15| as a basic model 
for the biological pattern formation as, for instance, tumor growth and bacterial colonies. In 
this model, new particles are sequentially added to the empty neighborhood of the cluster 



without overlap with previously aggregated particles [161, Il7|. Although the Eden model is 
unrealistic from the biological point of view, it produces compact aggregates with a nontrivial 
interface scaling usually analyzed through the interface width w ISj. Intensive numerical 
simulations indicate a power-law growth of the interface width with the time, w ~ t^, 
and exponent /3 = 1/3 jl4 . 17 , 19 . 20|. corresponding to the Kardar-Parisi-Zhang (KPZ) 
universality class |2l[ |. 

The DLA, BA, and Eden models can be implemented and simulated in a relatively easy 
way by constraining the particle positions to the sites of an underlying lattice. However, it 
is very well established that lattice anisotropy imposes strong effects on the cluster shape 



and scaling 22 



23 



24 



Althoug 



the anisotropy of on-lattice clusters 



1 some procedures have been proposed to remove 



26 



27 



28| . their successes were limited and off-lattice 



simulations impose themselves as a general framework for the investigation of the scaling 
properties and universality classes of these aggregates. Clearly, the aggregation of a large 
number of particles is necessary to reach the asymptotic behavior which, in turn, demands 
very efficient algorithms for large scale off-lattice simulations with rigorous statistical sam- 
pling. In this paper, we review several strategies used to optimize computer algorithms 
for off-lattice aggregates. Only those procedures which can be applied to general off-lattice 
simulations are focused here. More sophisticated but less general procedures, as conformal 



maps 



291], are avoided. Indeed, the conformal mapping is the most efficient strategy to 



simulate two-dimensional aggregates, but it cannot be used in higher dimensions. 



II. ALGORITHMS FOR OFF-LATTICE AGGREGATION 



In this section we present the description of distinct optimizations for two-dimensional 
clusters. The generalization for higher dimensions is straightforward. In all cases, simula- 
tions start with a single particle at the origin. 
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FIG. 1: (a) Schematic representation of the "optimized random trajectories", (b) A DLA aggregate 
and a mesh of cehs 2rint x 2rint- Long steps are forbidden in the gray boxes and allowed in the 
white ones. Also, two long steps are illustrated, (c) A zoom of the region inside the large square 
in (b). 



A. The trajectories 

Firstly, we describe optimizations for models in which particles of unitary diameter follow 
trajectories before stick to the aggregate, as is the case for the DLA and BA models. In 
both cases, the particles are released at random from a launching radius r/ larger than the 
cluster radius rmax and follow their trajectories up to touch the aggregate or cross a killing 
radius much larger than the system size. In the DLA, where particles follow discrete 
time random walks of unitary steps, a standard method is to allow the particles outside 
the launching circle take long random steps of length r^xt if these steps do not bring up 



a particle inside the launching circle, as illustrated in figure 1(a) An adequate choice is 
Text = max(r — r^ax — 1)? where r is the distance of the walker from the origin and a 
small tolerance 6 = 5 was used. Also, the Brownian walks in large empty areas in the inner 
region which delimits the cluster (r < r^ax) are very computer time consuming, specially for 



large aggregates. Ball and Brady [30| proposed a strategy which allows the particles inside 
the launching circle to take a long step of length Tint if they do not cross any part of the 
aggregate, as illustrated in figure 1(a), Similar procedures have been used in other works 



31 



32 



33|. 



In the BA model, the particles follow ballistic trajectories and the clusters do not exhibit 
large empty inner regions as in the DLA model. Hence, the trajectories can be efficiently 
implemented simply using a long step of size r^xt as in DLA model. An important difference 
between BA and DLA implementations is that in the first the launching radius should be as 
large as possible in order to avoid growth instabilities promoted by shadowing effects 34. l35| 
while in the DLA, this radius can be taken a few particle diameters larger than the cluster 
radius. 

A smart strategy to determine the length of the internal steps Vint is decisive for the 
algorithm efficiency. In order to accomplish this task, we define a square region of side 
L centered on the initial seed which delimits the entire aggregate. This region should be 
sufficiently large in order to guarantee that aggregate does not exceed its boundary. Then, 
the region is divided in a coarse-grained mesh with cells of size 2rj„t x 2rj„t as illustrated in 



figures 1(b) and 1(c) Each cell of the mesh is associated to an element oi a K x K square 
matrix A, where K = L/(2rj„t), which assumes 1 if the cell or one of its nearest or next- 
nearest neighbors contains any particle of the aggregate or assumes otherwise. The boxes 
depicted in gray [Aij = 1) are those in which the random walk can cross the cluster after a 
step of length Tint, since they contain or are adjacent to a part of the cluster. Consequently, 
long steps starting from gray boxes are forbidden. There are two options for a walker on a 
gray box: the particle executes a unitary step or tries a shorter step of length r^'^^, where 
1 < r^„j < Tint, using other auxiliary coarse-grained mesh A' with cells of size 2r-„^ x 2r^„j. 
Indeed, several auxiliary meshes can be used in order to maximize the efficiency. In this 
paper, we report simulation for 3 meshes with = 4, 8, and 16. 

The overlap between particles can occur after a unitary step if the preceding step brings 
the random walker at a distance from the cluster particle where it sticks lower than the 
unity. In this case, one just bring back the particle to the adjacent position along the 
opposite direction of the movement. 



B. Determination of the neighborhood 




(a) (b) 

FIG. 2: Illustration of the optimizations for off-lattice aggregation processes, (a) An auxiliary 
square lattice is used to determine when the walker is neighboring the cluster. The cluster particles 
are represented by black circles and their neighbors are depicted in gray, (b) A mesh with cells of 
size 4x4 used to restrict the search for contacts nearby the walker. 

The search mechanism for determining when and where the walker has contacted the 
aggregate represents the major time consuming step in large off-lattice simulations. The 
spatial coordinates of the particle belonging to the cluster are stored in one-dimensional 
arrays at the sequence of aggregation. So, the inspection of these arrays is performed 
whenever the walkers are in the nearby of the aggregate. If none optimization is adopted, 
all aggregated particles may be checked to verify if a contact occurred or not. At least 
three optimizations can be used. In the first and simplest one, we just verify the list in the 
reverse order in which the particles were added to the cluster, because the chance is larger 
for the aggregation to take place on the more external particles than on the inner ones. 
This procedure is considered default in this work. In the second one, particle positions are 
mapped on a square lattice by approximating their real coordinates to the nearest integer, 
producing an on-lattice cluster. In an auxiliary square lattice we label as occupied those 
sites belonging to the previous on-lattice cluster as well as their nearest and next-nearest 
neighbors. The search for contact is done only if the nearest integer coordinates of the 



walker represent an occupied site of the lattice Z. This procedure is schematically described 



in figure 2(a) In the third optimization procedure, a coarse-grained mesh W of cells with 
size £ X £ can be used to limit the verification to a region around the walker position. In 
this strategy, the cells are sequentially labeled by an index k = 1, 2, 3, • ■ • when they are 
occupied by a particle of the cluster for the first time. Also, the number of particles A^^ in 
the cells are stored. Finally, a third auxiliary one-dimensional array divided in blocks with 
i'^ elements is used to store the indexes of the particles in the arrays of coordinates. Each 
block is associated to a cell of the mesh. Once the analysis of the auxiliary square matrix 
Z have provided that the walker may be in contact with a particle of the cluster, the index 
k read in the mesh W is used to restrict the search for a contact in the array of coordinates 
using jr. The cell index of a walker at real coordinates {x,y) is given hj k = Wij, where 
i = nint(a;/^), j = nint (?//£), and nint(a;) function rounds x to the nearest integer. Indeed, 
the particles in the cell k are visited by varying the index of the array from n = + 1 
to n = Uk + Nk, where Uk = i'^ x {k — 1). Notice that the cell j of the mesh W and its 
neighbor cells should be verified to check the contacts on the cell edges. In the simulation 
results presented in the next section, £ = 4 was used. 



The Eden model 



The off-lattice simulation of the Eden model was proposed by Wang et al. 16| and 
improved by Ferreira and Alves 17|] as follows 



A particle with unitary diameter is chosen at random from a list of active ones (figure 



3(a)). A particle is considered active when a new one adjacent to it can be added to 
the aggregate without any overlap. 



Once an active particle was chosen (figure 3(b) ), its empty adjacent region, where there 
are no overlap between a new particle and those previously aggregated, is determined. 
A new particle is put in a direction randomly chosen among the allowed ones (figure 



3(c)). 



If the active particle does not have a growth region, it is labeled as inactive (figure 




(a) (b) (c) 

FIG. 3: Growth rules for the off-lattice Eden model. Active and inactive particles are represented 
by open and fullfiled discs, respectively, (a) A cluster and two active particles selected for the 
growth. The particle 1 has an empty region where a new adjacent particle can be added while 
the particle 2 does not. (b) The growth region adjacent to the particle 1 is shown as a dashed 
sector, (c) A new particle is added at a random direction in the growth region shown in (b) and 
the particle 2 is discarded from the list of active ones (both indicated by arrows) . 




(a) (b) (c) 



FIG. 4: (a) Eden cluster with 6000 particles. The border is represented by fullfiled symbols. Active 
particles for (b) standard and (c) optimized off-lattice algorithms for the Eden model are shown. 

In figure [3] the evolution rules are illustrated by two independent growth processes. Since 



the interest on Eden clusters is focused on the interface scaling, Ferreira and Alves [17 1 
introduced an optimization where any active cell inside a central core of radius is labeled 
as inactive. Since the inactivation of the particles near or belonging to the interface must be 
avoided, Vc = 0.8f was chosen, where f is the mean radius of the interface. This optimization 



was used only for r > 80a. In figure HI typical growth patterns with and without this last 



optimization, the corresponding borders 36|, and the active particles are illustrated. Finally, 
the optimizations described in sub-section III BI for determining the neighborhood of a particle 
can be used for the Eden model. 



III. SIMULATIONS 



All simulations were performed on the same computer, a Pentium IV 3.0 GHz with 2GB 
of RAM memory under Debian Linux operating system. One process was run by time. The 
algorithm codes were written in FORTRAN 90 language and compiled with the standard 
options of the Intel Fortran Compiler 9.1 37 1. 



A. Diffusion-limited aggregation 

Off-lattice DLA clusters with particles were grown using different combinations of the 
previously described optimizations. In all simulations, the launching and killing radius were 
taken as ri = Vmax + 5 and = lOOr^, respectively. In 1981, when Sander and Witten 
published their seminal work introducing the DLA model jc| without any optimization, the 
largest cluster generated on square lattices produced with computers of that age did not 
reach 4000 particles. Nowadays, this sort of simulation can be performed in a few minutes 
with any standard home computer. In tabled the CPU times spent in off-lattice simulations 
of a single cluster for some optimization schedules are listed. Also, CPU times are shown as 
functions of in figure [5l 

Simulations without optimizations become prohibitively long for relatively small aggre- 
gates. For example, a single cluster with 5 x 10^ particles demanded 10 days of simulations. 
If external steps are included in the original algorithm, for simplicity called by Oi, a great 
improvement of the efficiency is observed for very small clusters, but the simulations are 
also prohibitive for ~ 10^, since inner empty regions become of the same magnitude 
than the cluster size. Simulations become more than three orders of magnitude faster when 
the optimized neighborhood determination is included in Oi optimization, now called O2. 
Notice that the computational time grows approximately proportional to A^^ for both op- 
timizations Oq and O2. Simulations one order faster and CPU times growing slower with 
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TABLE I: Real CPU times in minutes for distinct optimizations applied to the DLA model. N is the 
number of particles; Oq refers to the algorithm with the default optimization where the backward 
inspetion of the coordinate arrays is used; Oi means that the long external steps of size r^xt were 
used; O2 means that external steps and optimized neighborhood were used simultaneously; O3 the 
previous optimizations plus the internal long steps of size rint (figures [1] and [2]) were adopted. The 
approximate dependence between CPU time and cluster size are indicated in the last line. 

increasing cluster size are performed when inner steps are included in O2 optimization. Also, 
notice that the computational time increases faster in Oi than in the others optimizations, 
but for large clusters Oq and Oi optimizations are expected to be equivalent due to the 
presence of large empty inner regions. 

B. Ballistic aggregation 

Off-lattice simulations of the BA model are very similar to the DLA model. The main 
difference is that the unitary steps performed by the walkers are in a fixed direction randomly 
chosen at the beginning of the ballistic walk. Also, the launching and killing radius used 
were r/ = lOOr^ax + 1000 and r/^ = ri + 10. In table [Tll the computational times for 
the same strategies used for DLA are listed. Like in the DLA model, long steps improve 




FIG. 5: CPU times as functions of the number of particles in the off-lattice DLA model for distinct 
optimization strategies. Lines are power fits. 

simulation efficiency for small clusters, but this gain decreases with increasing number of 
particles. However, optimized neighborhood determination provide a gain of three orders 
of magnitude. In figure [6] the CPU times are drawn as functions of A^. These times grow 
approximately as T ~ A^^ '', T ~ N^-^, and T ~ A^^ ° for Oq, Oi, and O2, respectively. 
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FIG. 6: CPU times as functions of the number of particles in the off- lattice BA model for distinct 
optimization strategies. Lines are power fits. 



N 


Oo 




O2 


1 X 


10^ 


2.81 X 10"^ 


1.51 X 10^2 


1.35 X 10-2 


2 X 


10^ 


6.45 X 10^1 


2.78 X 10^2 


2.09 X 10-2 


5 X 


10^ 


1.96 X 10° 


1.08 X 10-1 


2.27 X 10-2 


1 X 


10^ 


4.78 X 10° 


4.63 X 10-1 


3.62 X 10-2 


2 X 


10^ 


1.29 X 10^ 


2.03 X 10° 


6.10 X 10-2 


5 X 


10^ 


5.44 X 10^ 


1.43 X 10^ 


1.34 X 10-1 


1 X 


10^ 


1.67 X 10^ 


5.99 X 10^ 


2.52 X 10-1 


2 X 


10^ 


6.10 X 10^ 


2.85 X 10^ 


4.94 X 10-1 


5 X 


10^ 


3.50 X 10^ 


2.18 X 10^ 


1.22 X 10° 


1 X 


10^ 


— 


— 


2.43 X 10° 


2 X 


10^ 


— 


— 


5.01 X 10° 


5 X 


10^ 


— 


— 


1.29 X lOi 


5 X 


10^ 






2.61 X lOi 



TABLE II: Real CPU time in minutes for distinct optimizations applied to BA model. Optimiza- 
tions as in table HI 

C. Eden model 



The major challenge in off-lattice simulation of the Eden model is to determine which 
are the active cells. Since Eden model do not involve walkers, strategies as those of figures 
[Hand 



2 a^ 



do not have sense. But, an efficient determination of the empty neighborhood 
can be used as done for the DLA model. The original strategy proposed by Wang et al. 
[igI is called Eq and when the local search of neighbors is included, the model is denoted 
by El. CPU times are given in table UTTl and figure [71 The last algorithm overcomes the 
first one in three or more orders of magnitude. If a central core of particles is excluded from 
the list of active ones, the optimization E2, simulations becomes up to three times faster. 
Moreover, the efficiency gain increases with the number of particles. Indeed, CPU times 
grow approximately as T ~ A^2-5^ ^ A^i-^, and T ~ A^i'2 for Eo, Ei, and E2, respectively. 
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TABLE III: Eden Model Optimizations. Symbols £"0, Ei, and E2 described in text. 

lo'RTTq H 




N 

FIG. 7: CPU times as functions of the number of particles in the off-lattice Eden model using 
distinct optimization strategies. Lines are power fits. 

IV. SUMMARY 



Several optimizing strategies for the computer simulation of aggregation models dispersed 
throughout the literature were described in the present paper. It have been demonstrated 
that the combined implementation of such strategies can reduce in up to four order of 



magnitude the computer time demanded to perform large scale simulations of off-lattice 
aggregates with an increase of one order of magnitude in the allocated memory. Furthermore, 

these procedures can be applied to the simulations of other cluster growth processes beyond 
the traditional DLA, BA, and Eden models. 
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